
###########################
###########################
## Hungary data analysis
## Mares & Young (2017)
###########################
###########################

rm(list = ls())


## call packages - use install.packages() if necessary

library(lme4)
library(stargazer)
library(gtools)
library(data.table)
library(list)
library(xtable)
library(foreign)
library(car)
library(lmtest)
library(Amelia)
library(ggplot2)
library(arm)
library(nnet)


## set working directory 

setwd('/Users/leyou/Dropbox/Hungary_All/Analysis/core_rep_pkg')


## functions

st <- function(x) (x - mean(x, na.rm=T))/sd(x, na.rm=T)
gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}


## read in data

dat <- read.dta('data/data_final_20160411.dta')


##### 
## data prep
#####


## clean up age variables
dat$age_log <- log(dat$age)
dat$age_st <- (dat$age_log - mean(dat$age_log, na.rm=T))/sd(dat$age_log, na.rm=T)
dat$age_cat <- ifelse(dat$age < 30, 1, ifelse(dat$age < 50, 2, ifelse(dat$age < 65, 3, 4)))


## fix roma 
dat$roma <- ifelse(dat$roma==2, 0, dat$roma)


## standardize all variables
dat$core_st <- (dat$core - mean(dat$core, na.rm=T))/sd(dat$core, na.rm=T)
dat$age_cat_st <- (dat$age_cat - mean(dat$age_cat, na.rm=T))/sd(dat$age_cat, na.rm=T)
dat$female_st <- (dat$female - mean(dat$female, na.rm=T))/sd(dat$female, na.rm=T)
dat$roma_st <- (dat$roma - mean(dat$roma, na.rm=T))/sd(dat$roma, na.rm=T)
dat$hh_income_st <- (dat$hh_income - mean(dat$hh_income, na.rm=T))/sd(dat$hh_income, na.rm=T)
dat$empfull_st <- (dat$empfull - mean(dat$empfull, na.rm=T))/sd(dat$empfull, na.rm=T)
dat$close10_st <- (dat$close10 - mean(dat$close10, na.rm=T))/sd(dat$close10, na.rm=T)
dat$secure_st <- (dat$secure - mean(dat$secure, na.rm=T))/sd(dat$secure, na.rm=T)
dat$right_st <- (dat$right - mean(dat$right, na.rm=T))/sd(dat$right, na.rm=T)



## entitlements index
dat$entitles <- ifelse(dat$welfare_hh==1 | dat$debt==1 | dat$workfare==1, 1, 0)


## copartisanship 
dat$MayorParty2010 <- recode(dat$MayorParty2010, "'MZSP'='MSZP'")
dat$partym_voter <- ifelse(dat$party_id=="NA", NA, 
                            ifelse(dat$party_id==dat$MayorParty2010, 1, 0))
table(dat$partym_voter, dat$MayorParty2010)



#####
## Table 1: Summary stats
#####

## summary statistics table

vars = c('age', 'female', 'roma', 'hh_income', 'core', 'welfare_hh', 'debt', 'secret', 'fidesz_voter', 
         'right', 'voted')
levels = c(0.9, 0.95)
coefs = data.frame('mean_A' = rep(NA, length(vars)), 'mean_B' = rep(NA, length(vars)),
                   'diff' = rep(NA,length(vars)), 'p' = rep(NA, length(vars)), 'N' = rep(NA, length(vars)))
for (v in 1:length(vars)){
  for (j in 1:2) {
    c <- t.test(dat[,vars[v]] ~ dat[,'treat_welfare_pressure'],
                alternative = 'two.sided', conf.level = levels[j])
    coefs$diff[v] <- c$estimate[2] - c$estimate[1]
    coefs$mean_A[v] = c$estimate[1]; coefs$mean_B[v] = c$estimate[2]
    coefs$p[v] = c$p.value
    coefs$N[v] = length(na.omit(dat[,vars[v]]))
  }
}

rownames(coefs) = vars
xtable(coefs)


## other individual characteristics

prop.table(table(dat$hh_income, useNA = 'ifany'))
prop.table(table(dat$core, useNA = 'ifany'))
prop.table(table(dat$party_id[dat$party_id!="NA"]))
prop.table(table(dat$voted))
prop.table(table(dat$secret))


## community characteristics

table(dat$county)/length(dat$county)
table(dat$fidesz_mayor)/length(dat$fidesz_mayor)



#####
## Text p17: Test for design effects
#####

## design effects test

  ## call version of ict.test that allows NAs

source('functions/ict.test.R')
ict.test(y = dat$list_welfare_pressure, treat = dat$treat_welfare_pressure, J = 3)
ict.test(y = dat$list_lender_pressure, treat = dat$treat_lender_pressure, J = 3)
ict.test(y = dat$list_mayor_favor, treat = dat$treat_mayor_favor, J = 3)
ict.test(y = dat$list_vote_buying, treat = dat$treat_vote_buying, J = 3)

## control means and sd

stats <- data.frame('mean' = rep(NA, 4), 'sd' = rep(NA, 4))
stats$mean[1] <- mean(dat$list_welfare_pressure, na.rm=T)
stats$mean[2] <- mean(dat$list_lender_pressure, na.rm=T)
stats$mean[3] <- mean(dat$list_mayor_favor, na.rm=T)
stats$mean[4] <- mean(dat$list_vote_buying, na.rm=T)

stats$sd[1] <- sd(dat$list_welfare_pressure, na.rm=T)
stats$sd[2] <- sd(dat$list_lender_pressure, na.rm=T)
stats$sd[3] <- sd(dat$list_mayor_favor, na.rm=T)
stats$sd[4] <- sd(dat$list_vote_buying, na.rm=T)

stats


## percent control 0 or 3

table(dat$list_welfare_pressure[dat$treat_welfare_pressure==0]==0)/length(na.omit(dat$list_welfare_pressure[dat$treat_welfare_pressure==0]))
table(dat$list_lender_pressure[dat$treat_lender_pressure==0]==0)/length(na.omit(dat$list_lender_pressure[dat$treat_lender_pressure==0]))
table(dat$list_vote_buying[dat$treat_vote_buying==0]==0)/length(na.omit(dat$list_vote_buying[dat$treat_vote_buying==0]))
table(dat$list_mayor_favor[dat$treat_mayor_favor==0]==0)/length(na.omit(dat$list_mayor_favor[dat$treat_mayor_favor==0]))

table(dat$list_welfare_pressure[dat$treat_welfare_pressure==0]==3)/length(na.omit(dat$list_welfare_pressure[dat$treat_welfare_pressure==0]))
table(dat$list_lender_pressure[dat$treat_lender_pressure==0]==3)/length(na.omit(dat$list_lender_pressure[dat$treat_lender_pressure==0]))
table(dat$list_vote_buying[dat$treat_vote_buying==0]==3)/length(na.omit(dat$list_vote_buying[dat$treat_vote_buying==0]))
table(dat$list_mayor_favor[dat$treat_mayor_favor==0]==3)/length(na.omit(dat$list_mayor_favor[dat$treat_mayor_favor==0]))


## percent treatment 0 or 4

table(dat$list_welfare_pressure[dat$treat_welfare_pressure==1]==0)/length(na.omit(dat$list_welfare_pressure[dat$treat_welfare_pressure==1]))
table(dat$list_lender_pressure[dat$treat_lender_pressure==1]==0)/length(na.omit(dat$list_lender_pressure[dat$treat_lender_pressure==1]))
table(dat$list_vote_buying[dat$treat_vote_buying==1]==0)/length(na.omit(dat$list_vote_buying[dat$treat_vote_buying==1]))
table(dat$list_mayor_favor[dat$treat_mayor_favor==1]==0)/length(na.omit(dat$list_mayor_favor[dat$treat_mayor_favor==1]))

table(dat$list_welfare_pressure[dat$treat_welfare_pressure==1]==4)/length(na.omit(dat$list_welfare_pressure[dat$treat_welfare_pressure==1]))
table(dat$list_lender_pressure[dat$treat_lender_pressure==1]==4)/length(na.omit(dat$list_lender_pressure[dat$treat_lender_pressure==1]))
table(dat$list_vote_buying[dat$treat_vote_buying==1]==4)/length(na.omit(dat$list_vote_buying[dat$treat_vote_buying==1]))
table(dat$list_mayor_favor[dat$treat_mayor_favor==1]==4)/length(na.omit(dat$list_mayor_favor[dat$treat_mayor_favor==1]))




#####
## Figure 1: Prevalence of each strategy
#####

## prop of people overall who experience each strategy

  ## unweighted
vars = c('welfare_pressure', 'lender_pressure', 'mayor_favor', 'vote_buying')
levels = c(0.9, 0.95)
coefs = data.frame('coefs' = rep(NA,length(vars)), 'up90' = rep(NA,length(vars)), 'lo90' = rep(NA,length(vars)),
                   'up95' = rep(NA,length(vars)), 'lo95' = rep(NA,length(vars)))
for (v in 1:4){
  for (j in 1:2) {
    c <- t.test(dat[,paste0('list_',vars[v])] ~ dat[,paste0('treat_',vars[v])],
                alternative = 'two.sided', conf.level = levels[j])
    coefs$coefs[v] <- c$estimate[2] - c$estimate[1]
    if (j == 1){
      coefs$lo90[v] <- -c$conf.int[2]
      coefs$up90[v] <- -c$conf.int[1]
    }
    else {
      coefs$lo95[v] <- -c$conf.int[2]
      coefs$up95[v] <- -c$conf.int[1]
    }
  }
}


coefs = coefs*100
coefs$Strategy <- factor(c('Welfare
Pressure', 
'Lender
Pressure',
'Mayor
Favor',
'Vote
Buying'),
levels = c('Welfare
Pressure', 
           'Lender
Pressure',
           'Mayor
Favor',
           'Vote
Buying'))
coefs$Type <- c(rep("Negative", 2), rep('Positive', 2))


pdf('graphs/prop_strat.pdf')
ggplot(coefs) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = Strategy, ymin = lo95, ymax = up95, color = Type),
                lwd = 1, position = position_dodge(width = 1/2),
                width = 0) +
  geom_errorbar(aes(x = Strategy, ymin = lo90, ymax = up90, color = Type),
                lwd = 1, position = position_dodge(width = 1/2),
                width = 0.25) +
  geom_point(aes(x = Strategy, y = coefs, color = Type), 
             position = position_dodge(width=1/2), size = 4) +
  scale_y_continuous(name = "Percent of Respondents") + 
  theme_minimal() +
  theme(text = element_text(size = 20),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") 
dev.off()



#####
## Appendix Figure B2: Proportion by levels of closeness to party
#####

vars = c('welfare_pressure', 'lender_pressure', 'mayor_favor', 'vote_buying')
levels = c(0.9, 0.95)
coefs = data.frame('coefs' = rep(NA,4), 'up90' = rep(NA,4),
                   'lo90' = rep(NA,4),
                   'up95' = rep(NA,4), 'lo95' = rep(NA,4))

for (v in 1:length(vars)){
  for (i in 0:3) {
    for (j in 1:2) {
      c <- t.test(dat[,paste0('list_',vars[v])][dat$core==i] ~ dat[,paste0('treat_',vars[v])][dat$core==i],
                  alternative = 'two.sided', conf.level = levels[j])
      coefs$coefs[i+1] <- c$estimate[2] - c$estimate[1]
      if (j == 1){
        coefs$lo90[i+1] <- -c$conf.int[2]
        coefs$up90[i+1] <- -c$conf.int[1]
      }
      else {
        coefs$lo95[i+1] <- -c$conf.int[2]
        coefs$up95[i+1] <- -c$conf.int[1]
      }
    }
  }
  coefs$Core <- factor(c("No affiliation", "Not very 
close", "Somewhat
close", "Very close"),
                       levels = c("No affiliation", "Not very 
close", "Somewhat
close", "Very close"))
  assign(paste0('coefs',v), coefs)
}

pdf = c('welf', 'lend', 'may', 'buy')
list = list(coefs1, coefs2, coefs3, coefs4)
cols <- c(rep(gg_color_hue(2)[1], 2), rep(gg_color_hue(2)[2], 2))

for (i in 1:length(vars)){
  coefs <- list[[i]]
  g <- ggplot(coefs) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_errorbar(aes(x = Core, ymin = lo95, ymax = up95),
                  lwd = 1, position = position_dodge(width = 1/2),
                  width = 0, color = cols[i]) +
    geom_errorbar(aes(x = Core, ymin = lo90, ymax = up90),
                  lwd = 1, position = position_dodge(width = 1/2),
                  width = 0.25, color = cols[i]) +
    geom_point(aes(x = Core, y = coefs), 
               position = position_dodge(width=1/2), size = 4, color = cols[i]) +
    scale_y_continuous(name = "Percent of Respondents") + 
    theme_minimal() +
    theme(text = element_text(size = 20),
          axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = "none") 
  pdf(paste0('graphs/prop_core_', pdf[i], '.pdf'))
  print(g)
  dev.off()
}



#####
## Table 2: Closeness to party and clientelism 
#####

spec1 <- 'core_st'
spec2 <- 'core_st + age_cat_st + roma + female + hh_income_st'

dat$treat <- dat$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
             data = dat, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
             data = dat, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
             data = dat, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
             data = dat, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')


## Output table 2 

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(p5$par.treat)*4),
                        'mod2' = rep(NA, length(p5$par.treat)*4),
                        'mod3' = rep(NA, length(p5$par.treat)*4),
                        'mod4' = rep(NA, length(p5$par.treat)*4),
                        'mod5' = rep(NA, length(p5$par.treat)*4),
                        'mod6' = rep(NA, length(p5$par.treat)*4),
                        'mod7' = rep(NA, length(p5$par.treat)*4),
                        'mod8' = rep(NA, length(p5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(p5$par.treat)*4,2)] <- c(paste0(names(p4$par.treat),".t"), 
                                paste0(names(p4$par.control),".c"))

ptmax <- length(p5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))



#####
## Table 3: Closeness to party and access to entitlements
#####

d1 <- lm(debt ~ core_st, data = dat)
d2 <- lm(debt ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)
w1 <- lm(welfare_hh ~ core_st, data = dat)
w2 <- lm(welfare_hh ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)
l1 <- lm(workfare ~ core_st, data = dat)
l2 <- lm(workfare ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)
e1 <- lm(entitles ~ core_st, data = dat)
e2 <- lm(entitles ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)

stargazer(d1, d2, w1, w2, l1, l2, e1, e2,
          no.space = T, 
          covariate.labels = c('Core', 'Age', 'Roma', 'Female', 'Income', '(Intercept)'), 
          keep.stat = c('n', 'rsq'))



#####
## Table 4: Access to entitlements and election-time clientelism
#####

spec1 <- 'entitles'
spec2 <- 'entitles + age_cat_st + roma + female + hh_income_st'


dat$treat <- dat$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')


## table like stargazer

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(p5$par.treat)*4),
                        'mod2' = rep(NA, length(p5$par.treat)*4),
                        'mod3' = rep(NA, length(p5$par.treat)*4),
                        'mod4' = rep(NA, length(p5$par.treat)*4),
                        'mod5' = rep(NA, length(p5$par.treat)*4),
                        'mod6' = rep(NA, length(p5$par.treat)*4),
                        'mod7' = rep(NA, length(p5$par.treat)*4),
                        'mod8' = rep(NA, length(p5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(p5$par.treat)*4,2)] <- c(paste0(names(p4$par.treat),".t"), 
                                                    paste0(names(p4$par.control),".c"))

ptmax <- length(p5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))



#####
## Table 5: are copartisans of the mayor targeted?
#####

table(dat$party_id[dat$partym_voter==1])


spec1 <- 'partym_voter'
spec2 <- 'partym_voter + age_cat_st + roma + female + hh_income_st'

dat$treat <- dat$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')


  ## table like stargazer

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(p5$par.treat)*4),
                        'mod2' = rep(NA, length(p5$par.treat)*4),
                        'mod3' = rep(NA, length(p5$par.treat)*4),
                        'mod4' = rep(NA, length(p5$par.treat)*4),
                        'mod5' = rep(NA, length(p5$par.treat)*4),
                        'mod6' = rep(NA, length(p5$par.treat)*4),
                        'mod7' = rep(NA, length(p5$par.treat)*4),
                        'mod8' = rep(NA, length(p5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(p5$par.treat)*4,2)] <- c(paste0(names(p4$par.treat),".t"), 
                                                    paste0(names(p4$par.control),".c"))

ptmax <- length(p5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))



#####
## Appendix B.1: Robustness to categorical version of core
#####

dat$core_cat <- as.factor(dat$core)
dat$hh_income_cat <- as.factor(dat$hh_income)


## Table B1 - core and election-time clientelism

spec1 <- 'core_cat'
spec2 <- 'core_cat + age_cat_st + roma + female + hh_income_cat'

dat$treat <- dat$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')


## table like stargazer

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(p5$par.treat)*4),
                        'mod2' = rep(NA, length(p5$par.treat)*4),
                        'mod3' = rep(NA, length(p5$par.treat)*4),
                        'mod4' = rep(NA, length(p5$par.treat)*4),
                        'mod5' = rep(NA, length(p5$par.treat)*4),
                        'mod6' = rep(NA, length(p5$par.treat)*4),
                        'mod7' = rep(NA, length(p5$par.treat)*4),
                        'mod8' = rep(NA, length(p5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(p5$par.treat)*4,2)] <- c(paste0(names(p4$par.treat),".t"),
                                                    paste0(names(p4$par.control),".c"))

ptmax <- length(p5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))


## Figure B2 - graph coefficients in categorical specs w controls

models <- list(p4, p5, n4, n5)
coefs <- data.frame('Core' = factor(c('Not very
close', 'Somewhat
close', 'Very
close'), levels = c('Not very
close', 'Somewhat
close', 'Very
close')))

for (i in 1:length(models)){
  mod <- models[[i]]
  coefs$coefs <- c(mod$par.treat[grep('core_cat[0-9]',names(mod$par.treat))])
  coefs$ses <- c(mod$se.treat[grep('core_cat[0-9]',names(mod$se.treat))])
  coefs$up <- coefs$coefs + 1.96 * coefs$ses; lo <- coefs$coefs - 1.96 * coefs$ses
  coefs$up90 <- coefs$coefs + 1.64 * coefs$ses; lo90 <- coefs$coefs - 1.64 * coefs$ses
  colors <- if (i > 2){gg_color_hue(2)[1]} else {{gg_color_hue(2)[2]}}
  g <- ggplot(coefs) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_errorbar(aes(x = Core, ymin = lo, ymax = up),
                  lwd = 1, position = position_dodge(width = 1/2),
                  width = 0, color = colors) +
    geom_errorbar(aes(x = Core, ymin = lo90, ymax = up90),
                  lwd = 1, position = position_dodge(width = 1/2),
                  width = 0.25, color = colors) +
    geom_point(aes(x = Core, y = coefs), 
               position = position_dodge(width=1/2), size = 4, color = colors) +
    scale_y_continuous(name = "Coefficient") + 
    theme_minimal() +
    theme(text = element_text(size = 20),
          axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = "none") 
  pdf(paste0("graphs/coefs_cat_", i, ".pdf"))
  print(g)
  dev.off()
}



## Table B2: Closeness to party and access to entitlements - categorical

d1 <- lm(debt ~ core_cat, data = dat)
d2 <- lm(debt ~ core_cat + age_cat_st + roma + female + hh_income_st, data = dat)
w1 <- lm(welfare_hh ~ core_cat, data = dat)
w2 <- lm(welfare_hh ~ core_cat + age_cat_st + roma + female + hh_income_st, data = dat)
l1 <- lm(workfare ~ core_cat, data = dat)
l2 <- lm(workfare ~ core_cat + age_cat_st + roma + female + hh_income_st, data = dat)
e1 <- lm(entitles ~ core_cat, data = dat)
e2 <- lm(entitles ~ core_cat + age_cat_st + roma + female + hh_income_st, data = dat)

stargazer(d1, d2, w1, w2, l1, l2, e1, e2,
          no.space = T, 
          covariate.labels = c('Core - Cat 1', 'Core - Cat 2', 'Core - Cat 3',
                               'Age', 'Roma', 'Female', 'Income', '(Intercept)'), 
          keep.stat = c('n', 'rsq'))



## Figure B3: graph coefficients in categorical

models <- list(d2, w2, l2, e2)
coefs <- data.frame('Core' = factor(c('Not very
close', 'Somewhat
close', 'Very
close'), levels = c('Not very
close', 'Somewhat
close', 'Very
close')))

for (i in 1:length(models)){
  mod <- summary(models[[i]])$coef
  coefs$coefs <- c(mod[grep('core_cat[0-9]',rownames(mod)), 1])
  coefs$ses <- c(mod[grep('core_cat[0-9]',rownames(mod)), 2])
  coefs$up <- coefs$coefs + 1.96 * coefs$ses; lo <- coefs$coefs - 1.96 * coefs$ses
  coefs$up90 <- coefs$coefs + 1.64 * coefs$ses; lo90 <- coefs$coefs - 1.64 * coefs$ses
  colors <- gg_color_hue(2)[1]
  g <- ggplot(coefs) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_errorbar(aes(x = Core, ymin = lo, ymax = up),
                  lwd = 1, position = position_dodge(width = 1/2),
                  width = 0, color = colors) +
    geom_errorbar(aes(x = Core, ymin = lo90, ymax = up90),
                  lwd = 1, position = position_dodge(width = 1/2),
                  width = 0.25, color = colors) +
    geom_point(aes(x = Core, y = coefs), 
               position = position_dodge(width=1/2), size = 5, color = colors) +
    scale_y_continuous(name = "Coefficient") + 
    theme_minimal() +
    theme(text = element_text(size = 22),
          axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = "none") 
  pdf(paste0("graphs/coefs_entitles_cat_", i, ".pdf"))
  print(g)
  dev.off()
}




#####
## Appendix B.2, B.3: Robustness to alternative measures of core
#####

## Figure B4: party id and ideology

coefs <- c(mean(dat$right[dat$fidesz_voter==1],na.rm=T),
           mean(dat$right[dat$jobbik_voter==1],na.rm=T),
           mean(dat$right[dat$soc_voter==1],na.rm=T))
ses <- c(sd(dat$right[dat$fidesz_voter==1],na.rm=T)/sqrt(length(dat$right[dat$fidesz_voter==1])),
         sd(dat$right[dat$jobbik_voter==1],na.rm=T)/sqrt(length(dat$right[dat$jobbik_voter==1])),
         sd(dat$right[dat$soc_voter==1],na.rm=T)/sqrt(length(dat$right[dat$soc_voter==1])))
n <- c(length(na.omit(dat$right[dat$fidesz_voter==1])),
       length(na.omit(dat$right[dat$jobbik_voter==1])),
       length(na.omit(dat$right[dat$soc_voter==1])))
colors = c('red', 'green', 'blue')

pdf("graphs/ideology_party.pdf")
plot(y = seq(1,3,1), x = coefs, col = colors, pch = 16,
     xlim = c(-.5,.15), 
     ylim = c(-5,15),
     xlab = 'Ideology on Left-Right Scale', 
     ylab = "", yaxt = "n")
for (i in 1:3){
  lines(y = rep(i, 2), x = c(coefs[i]-1.96*ses[i], coefs[i]+1.96*ses[i]), 
        col = colors[i], lwd = 2)
  lines(y = rep(i, 2), x = c(coefs[i]-1.64*ses[i], coefs[i]+1.64*ses[i]), 
        col = colors[i], lwd = 4)
}
legend(col = c('blue', 'red', 'green'), pch = 16, ncol = 3,
       legend = c('Socialist', 'Fidesz', 'Jobbik'),
       'top')
text(x = coefs, y = seq(1,3,1)-2, paste0('N=',n))
dev.off()


## Table B3: Ideology and closeness to party

dat$party_voter <- as.factor(ifelse(dat$fidesz_voter==1, 'fidesz',
                                    ifelse(dat$jobbik_voter==1, 'jobbik',
                                           ifelse(dat$soc_voter==1, 'mszp',
                                                  ifelse(dat$noparty_voter==1, 'none', 'other')))))

dat$dats <- ifelse((dat$soc_voter==1 | dat$core==0) & is.na(dat$right_st)==F, 1, 0)
dat$dats <- ifelse(is.na(dat$dats), 0, dat$dats)
dat$datf <- ifelse((dat$fidesz_voter==1 | dat$core==0) & is.na(dat$right_st)==F, 1, 0)
dat$datf <- ifelse(is.na(dat$datf), 0, dat$datf)
dat$datj <- ifelse((dat$jobbik_voter==1 | dat$core==0) & is.na(dat$right_st)==F, 1, 0)
dat$datj <- ifelse(is.na(dat$datj), 0, dat$datj)

i1a <- lm(core_st ~ right_st,
          data = dat[dat$dats==1,])
i1b <- lm(core_st ~ right_st + female + hh_income_st + roma + age_st,
          data = dat[dat$dats==1,])

i2a <- lm(core_st ~ right_st,
          data = dat[dat$datf==1,])
i2b <- lm(core_st ~ right_st + female + hh_income_st + roma + age_st,
          data = dat[dat$datf==1,])

i3a <- lm(core_st ~ right_st, 
          data = dat[dat$datj==1,])
i3b <- lm(core_st ~ right_st + female + hh_income_st + roma + age_st, 
          data = dat[dat$datj==1,])

stargazer(i1a, i1b, i2a, i2b, i3a, i3b, 
          covariate.labels = c('Right', 'Female', 'Income', 'Roma', 'Age', '(Intercept)'),
          no.space = T, keep.stat = c('n', 'rsq'))


## Table B4: Alternative measures of closeness to party and clientelism targeting

  ## create predicted measure of core by ideology and party

i1 <- lm(core_st ~ right_st + party_voter + right_st:party_voter, data = dat)
inc <- complete.cases(dat[,colnames(i1$model)])
dat$core_pred2[inc==T] <- predict(i1)
dat$core_pred2_st <- st(dat$core_pred2)

    ## check how well the predicted value corresponds to the actual level of core

summary(lm(core_pred2_st ~ core_st, data = dat))
summary(lm(core_pred2_st ~ core_st, data = dat[dat$datf==1,]))
summary(lm(core_pred2_st ~ core_st, data = dat[dat$datj==1,]))
summary(lm(core_pred2_st ~ core_st, data = dat[dat$dats==1,]))


  ## create binary measure of core

dat$core_bin <- ifelse(dat$core > 0, 1, 0)



  ## specifications

spec1 <- 'core_temp'
spec2 <- 'core_temp + age_cat_st + roma + female + hh_income_st'

dat$core_temp <- dat$core_st

dat$treat <- dat$treat_vote_buying
p1a <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p1b <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p1c <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p1d <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n1a <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n1b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n1c <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n1d <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')

dat$core_temp <- dat$core_bin

dat$treat <- dat$treat_vote_buying
p2a <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p2b <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p2c <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p2d <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n2a <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n2b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n2c <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n2d <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')

dat$core_temp <- dat$core_pred2_st

dat$treat <- dat$treat_vote_buying
p3a <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p3b <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p3c <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p3d <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n3a <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n3b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n3c <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n3d <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
              data = dat, treat = "treat", J = 3, method = 'lm')


  ## table like stargazer

models1 <- list(p1a, p1b, p1c, p1d, n1a, n1b, n1c, n1d)
models2 <- list(p2a, p2b, p2c, p2d, n2a, n2b, n2c, n2d)
models3 <- list(p3a, p3b, p3c, p3d, n3a, n3b, n3c, n3d)
models <- list(models1, models2, models3)
tab <- cbind.data.frame('mod1' = rep(NA, 6),
                        'mod2' = rep(NA, 6),
                        'mod3' = rep(NA, 6),
                        'mod4' = rep(NA, 6),
                        'mod5' = rep(NA, 6),
                        'mod6' = rep(NA, 6),
                        'mod7' = rep(NA, 6),
                        'mod8' = rep(NA, 6))
modinf <- rbind.data.frame('N' = rep(NA, length(models)),
                           'N' = rep(NA, length(models)),
                           'N' = rep(NA, length(models)))
rownames(tab)[c(1,3,5)] <- c('Core', 'Core (Binary)', 'Core (Predicted)')

for (j in 1:3){
  model <- models[[j]]
  for (i in 1:8){
    mod <- model[[i]]
    
    mod$pval <- coeftest(mod)[,4]
    mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
    tab[j*2-1, i] <- paste0(round(mod$par.treat['core_temp'], 3), mod$stars['sensitive.core_temp'])
    tab[j*2, i] <- paste0("(",round(mod$se.treat['core_temp'],2),")")
    modinf[j, i] <- dim(mod$data)[1]
  }
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))



## Table B5: Entitlements and robustness to alternative core measures

  ## original core

d1 <- lm(debt ~ core_st, data = dat)
d2 <- lm(debt ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)
w1 <- lm(welfare_hh ~ core_st, data = dat)
w2 <- lm(welfare_hh ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)
l1 <- lm(workfare ~ core_st, data = dat)
l2 <- lm(workfare ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)
e1 <- lm(entitles ~ core_st, data = dat)
e2 <- lm(entitles ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat)

stargazer(d1, d2, w1, w2, l1, l2, e1, e2,
          no.space = T, 
          covariate.labels = c('Core', 'Age', 'Roma', 'Female', 'Income', '(Intercept)'),
          keep.stat = c('n', 'rsq'))


  ## binary core

d1 <- lm(debt ~ core_bin, data = dat)
d2 <- lm(debt ~ core_bin + age_cat_st + roma + female + hh_income_st, data = dat)
w1 <- lm(welfare_hh ~ core_bin, data = dat)
w2 <- lm(welfare_hh ~ core_bin + age_cat_st + roma + female + hh_income_st, data = dat)
l1 <- lm(workfare ~ core_bin, data = dat)
l2 <- lm(workfare ~ core_bin + age_cat_st + roma + female + hh_income_st, data = dat)
e1 <- lm(entitles ~ core_bin, data = dat)
e2 <- lm(entitles ~ core_bin + age_cat_st + roma + female + hh_income_st, data = dat)

stargazer(d1, d2, w1, w2, l1, l2, e1, e2,
          no.space = T, 
          covariate.labels = c('Core', 'Age', 'Roma', 'Female', 'Income', '(Intercept)'),
          keep.stat = c('n', 'rsq'))


  ## predicted core by ideology

d1 <- lm(debt ~ core_pred2_st, data = dat)
d2 <- lm(debt ~ core_pred2_st + age_cat_st + roma + female + hh_income_st, data = dat)
w1 <- lm(welfare_hh ~ core_pred2_st, data = dat)
w2 <- lm(welfare_hh ~ core_pred2_st + age_cat_st + roma + female + hh_income_st, data = dat)
l1 <- lm(workfare ~ core_pred2_st, data = dat)
l2 <- lm(workfare ~ core_pred2_st + age_cat_st + roma + female + hh_income_st, data = dat)
e1 <- lm(entitles ~ core_pred2_st, data = dat)
e2 <- lm(entitles ~ core_pred2_st + age_cat_st + roma + female + hh_income_st, data = dat)

stargazer(d1, d2, w1, w2, l1, l2, e1, e2,
          no.space = T, 
          covariate.labels = c('Core', 'Age', 'Roma', 'Female', 'Income', '(Intercept)'),
          keep.stat = c('n', 'rsq'))





#####
## Appendix C.1: less or more than they deserve?
#####

table(dat$empfull, dat$workfare)
table(dat$empfull, dat$welfare_hh)

dat$age_ret <- ifelse(dat$age > 62, 1, 0)

table(dat$age_ret, dat$workfare)
table(dat$age_ret, dat$welfare_hh)

prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$core==0])) # 23% of non-affiliated voters get welfare
prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$core==1])) # 26% get welfare when core=1
prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$core==2])) # 29% get welfare when core=2
prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$core==3])) # 30% get welfare when core=3

  ## proportion of workfare eligible voters that get welfare
t1 <- data.frame('no' = rep(NA, 5), 'yes' = rep(NA, 5))
t1[1,] <- prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$age <= 62])) 
t1[2,] <- prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$age <= 62 & dat$core==0])) 
t1[3,] <- prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$age <= 62  & dat$core==1])) 
t1[4,] <- prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$age <= 62  & dat$core==2])) 
t1[5,] <- prop.table(table(dat$welfare_hh[dat$empfull==0 & dat$age <= 62  & dat$core==3])) 
t1$x <- c('all', 'No affil.', 'Not very', 'Somewhat', 'Very')

  ## proportion of workfare eligible voters that get workfare
t2 <- data.frame('no' = rep(NA, 5), 'yes' = rep(NA, 5))
t2[1,] <- prop.table(table(dat$workfare[dat$age <= 62])) 
t2[2,] <- prop.table(table(dat$workfare[dat$age <= 62 & dat$core==0])) 
t2[3,] <- prop.table(table(dat$workfare[dat$age <= 62  & dat$core==1])) 
t2[4,] <- prop.table(table(dat$workfare[dat$age <= 62  & dat$core==2])) 
t2[5,] <- prop.table(table(dat$workfare[dat$age <= 62  & dat$core==3])) 
t2$x <- c('all', 'No affil.', 'Not very', 'Somewhat', 'Very')

  ## proportion of workfare ineligible voters that get welfare
prop.table(table(dat$welfare_hh[(dat$empfull==1 | dat$age > 62)]))
prop.table(table(dat$welfare_hh[(dat$empfull==1 | dat$age > 62) & dat$core==0])) 
prop.table(table(dat$welfare_hh[(dat$empfull==1 | dat$age > 62)  & dat$core==1]))
prop.table(table(dat$welfare_hh[(dat$empfull==1 | dat$age > 62)  & dat$core==2])) 
prop.table(table(dat$welfare_hh[(dat$empfull==1 | dat$age > 62)  & dat$core==3])) 

  ## proportion of elderly that get workfare
table(dat$workfare, dat$age > 62)

  ## proportion not fully employed
prop.table(table(dat$empfull))
prop.table(table(dat$empfull[dat$age <= 62]))


## Figure C5: plot proportion eligible

pdf("graphs/prop_elig_welf.pdf")
ggplot(t1[2:5,], aes(x, yes)) +
  geom_bar(stat = 'identity', fill = c(gg_color_hue(2)[1])) +
  theme_minimal() +
  theme(text = element_text(size=30), 
        axis.text.x = element_text(angle=45, vjust=0.95, hjust=1),
        legend.position = "none") +
  scale_y_continuous(labels = scales::percent, limit = c(0,0.6)) +
  ylab("Proportion") + xlab("Closeness to Party") 
dev.off()

pdf("graphs/prop_elig_work.pdf")
ggplot(t2[2:5,], aes(x, yes)) +
  geom_bar(stat = 'identity', fill = c(gg_color_hue(2)[1])) +
  theme_minimal() +
  theme(text = element_text(size=30), 
        axis.text.x = element_text(angle=45, vjust=0.95, hjust=1),
        legend.position = "none") +
  scale_y_continuous(labels = scales::percent, limit = c(0,0.2)) +
  ylab("Proportion") + xlab("Closeness to Party") 
dev.off()



#####
## Appendix C.2: Disaggregated entitlements results
#####

## Table C6: Who gets entitlements? Individual measures

spec1 <- 'welfare_hh + workfare + debt'
spec2 <- 'welfare_hh + workfare + debt + age_cat_st + roma + female + hh_income_st'

dat$treat <- dat$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')

dat$treat <- dat$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat, treat = "treat", J = 3, method = 'lm')


## table like stargazer

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(n5$par.treat)*4),
                        'mod2' = rep(NA, length(n5$par.treat)*4),
                        'mod3' = rep(NA, length(n5$par.treat)*4),
                        'mod4' = rep(NA, length(n5$par.treat)*4),
                        'mod5' = rep(NA, length(n5$par.treat)*4),
                        'mod6' = rep(NA, length(n5$par.treat)*4),
                        'mod7' = rep(NA, length(n5$par.treat)*4),
                        'mod8' = rep(NA, length(n5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(n5$par.treat)*4,2)] <- c(paste0(names(n5$par.treat),".t"), 
                                                    paste0(names(n5$par.control),".c"))

ptmax <- length(n5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))




#####
## Appendix D: Copartisanship
#####

## do copartisans gets entitlements?

d1 <- lm(debt ~ partym_voter, data = dat)
d2 <- lm(debt ~ partym_voter + age_cat_st + roma + female + hh_income_st, data = dat)
w1 <- lm(welfare_hh ~ partym_voter, data = dat)
w2 <- lm(welfare_hh ~ partym_voter + age_cat_st + roma + female + hh_income_st, data = dat)
l1 <- lm(workfare ~ partym_voter, data = dat)
l2 <- lm(workfare ~ partym_voter + age_cat_st + roma + female + hh_income_st, data = dat)
e1 <- lm(entitles ~ partym_voter, data = dat)
e2 <- lm(entitles ~ partym_voter + age_cat_st + roma + female + hh_income_st, data = dat)

stargazer(d1, d2, w1, w2, l1, l2, e1, e2,
          no.space = T, 
          covariate.labels = c('Mayor Copartisan', 'Age', 'Roma', 'Female', 'Income', '(Intercept)'),
          keep.stat = c('n', 'rsq'))



#####
## Appendix E: Missingness
#####

## Table E8: Correlates of missingness

dat$core_miss <- ifelse(is.na(dat$core), 1, 0)
dat$hh_income_miss <- ifelse(is.na(dat$hh_income), 1, 0)

m1a <- lm(core_miss ~ age_cat_st + roma + female + right_st, data = dat)
m1b <- lm(core_miss ~ age_cat_st + roma + female + right_st + debt + welfare_hh, data = dat)
m2a <- lm(hh_income_miss ~ age_cat_st + roma + female + right_st, data = dat)
m2b <- lm(hh_income_miss ~ age_cat_st + roma + female + right_st + debt + welfare_hh, data = dat)

stargazer(m1a, m1b, m2a, m2b, 
          no.space = T,
          covariate.labels = c('Age', 'Roma', 'Female', 'Conservative', 'Debt', 'Welfare', '(Intercept)'),
          keep.stat = c('n', 'rsq'))



## imputation

set.seed(84774)
a <- subset(dat, select =c(hh_income, d1, d12, d13, d14, debt, welfare_hh, empfull, workfare, party_id,
                           core, op_welfare_roma, op_crime, op_poor, op_school, 
                           voted, secret, female, age, roma, municipalityname,
                           treat_welfare_pressure, treat_mayor_favor, treat_lender_pressure, treat_vote_buying, 
                           list_mayor_favor, list_vote_buying, list_welfare_pressure, list_lender_pressure,
                           MayorParty2010))
a$party_id[a$party_id=="NA"] <- NA
a$party_id <- as.factor(a$party_id)
a$municipalityname <- as.factor(a$municipalityname)
a$d1[a$d1=="Refused"] <- NA
a$d1 <- as.factor(a$d1)
setnames(a, 'd1', 'emp_status')
a$d12 <- gsub(" *$", "", as.character(a$d12))
a$d12[a$d12 %in% c("Don’t know", "Refused")] <- NA
setnames(a, 'd12', 'feels_discrim')
a$feels_discrim <- as.factor(a$feels_discrim)
a$d13 <- gsub(" *$", "", as.character(a$d13))
a$d13[a$d13 %in% c("Don’t know", "Refused")] <- NA
setnames(a, 'd13', 'lend_500')
a$lend_500 <- as.factor(a$lend_500)
a$d14 <- gsub(" *$", "", as.character(a$d14))
a$d14[a$d14 %in% c("Don’t know", "Refused")] <- NA
setnames(a, 'd14', 'recip_500')
a$recip_500 <- as.factor(a$recip_500)
a$workfare=NULL; a$empfull=NULL

a.out <- amelia(a, ords = c('party_id', 'municipalityname', 'debt', 'emp_status',
                            'feels_discrim', 'lend_500', 'recip_500',
                            'core', 'hh_income', 'welfare_hh', 
                            'op_welfare_roma', 'op_crime', 
                            'op_poor', 'op_school', 'voted', 'secret', 
                            'female', 'roma', 
                            'list_mayor_favor', 'list_vote_buying', 
                            'list_welfare_pressure', 'list_lender_pressure'),
                idvars = c('treat_welfare_pressure', 'treat_lender_pressure',
                           'treat_vote_buying', 'treat_mayor_favor',
                           'MayorParty2010'))


## prep imputed data

dat.imp <- a.out$imputations[[1]]
dat.imp$core_st <- st(dat.imp$core)
dat.imp$age_cat <- ifelse(dat.imp$age < 30, 1, ifelse(dat.imp$age < 50, 2, ifelse(dat.imp$age < 65, 3, 4)))
dat.imp$age_cat_st <- st(dat.imp$age_cat)
dat.imp$hh_income_st <- st(dat.imp$hh_income)
dat.imp$MayorParty2010 <- as.character(recode(dat.imp$MayorParty2010, "'MZSP'='MSZP'"))
dat.imp$party_id <- as.character(dat.imp$party_id)
dat.imp$partym_voter <- ifelse(dat.imp$party_id=="NA", NA, 
                            ifelse(dat.imp$party_id==dat.imp$MayorParty2010, 1, 0))
dat.imp$workfare <- ifelse(dat.imp$emp_status=="Work in a welfare work program", 1, 0)
table(dat.imp$partym_voter, dat.imp$party_id)
table(dat.imp$MayorParty2010, dat.imp$party_id)


## Table E9 - imputed core and election-time clientelism

spec1 <- 'core_st'
spec2 <- 'core_st + age_cat_st + roma + female + hh_income_st'

dat.imp$treat <- dat.imp$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')

dat.imp$treat <- dat.imp$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')


  ## stargazer table

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(p5$par.treat)*4),
                        'mod2' = rep(NA, length(p5$par.treat)*4),
                        'mod3' = rep(NA, length(p5$par.treat)*4),
                        'mod4' = rep(NA, length(p5$par.treat)*4),
                        'mod5' = rep(NA, length(p5$par.treat)*4),
                        'mod6' = rep(NA, length(p5$par.treat)*4),
                        'mod7' = rep(NA, length(p5$par.treat)*4),
                        'mod8' = rep(NA, length(p5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(p5$par.treat)*4,2)] <- c(paste0(names(p4$par.treat),".t"), 
                                                    paste0(names(p4$par.control),".c"))

ptmax <- length(p5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))


## Table E10: imputed core and access to entitlements

dat.imp$entitles <- ifelse(dat.imp$welfare_hh==1 | dat.imp$debt==1 | dat.imp$workfare==1, 1, 0)

d1 <- lm(debt ~ core_st, data = dat.imp)
d2 <- lm(debt ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat.imp)
w1 <- lm(welfare_hh ~ core_st, data = dat.imp)
w2 <- lm(welfare_hh ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat.imp)
l1 <- lm(workfare ~ core_st, data = dat.imp)
l2 <- lm(workfare ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat.imp)
e1 <- lm(entitles ~ core_st, data = dat.imp)
e2 <- lm(entitles ~ core_st + age_cat_st + roma + female + hh_income_st, data = dat.imp)

stargazer(d1, d2, w1, w2, l1, l2, e1, e2,
          no.space = T, 
          covariate.labels = c('Core', 'Age', 'Roma', 'Female', 'Income', '(Intercept)'), 
          keep.stat = c('n', 'rsq'))



## Table E11: are those who get entitlements more targeted?

spec1 <- 'entitles'
spec2 <- 'entitles + age_cat_st + roma + female + hh_income_st'

dat.imp$treat <- dat.imp$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')

dat.imp$treat <- dat.imp$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')


  ## table like stargazer

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(p5$par.treat)*4),
                        'mod2' = rep(NA, length(p5$par.treat)*4),
                        'mod3' = rep(NA, length(p5$par.treat)*4),
                        'mod4' = rep(NA, length(p5$par.treat)*4),
                        'mod5' = rep(NA, length(p5$par.treat)*4),
                        'mod6' = rep(NA, length(p5$par.treat)*4),
                        'mod7' = rep(NA, length(p5$par.treat)*4),
                        'mod8' = rep(NA, length(p5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(p5$par.treat)*4,2)] <- c(paste0(names(p4$par.treat),".t"), 
                                                    paste0(names(p4$par.control),".c"))

ptmax <- length(p5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))



## Table E12: copartisanship with the mayor 

spec1 <- 'partym_voter'
spec2 <- 'partym_voter + age_cat_st + roma + female + hh_income_st'

dat.imp$treat <- dat.imp$treat_vote_buying
p4b <- ictreg(formula(paste0('list_vote_buying ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
p4 <- ictreg(formula(paste0('list_vote_buying ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')
p5b <- ictreg(formula(paste0('list_mayor_favor ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
p5 <- ictreg(formula(paste0('list_mayor_favor ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')

dat.imp$treat <- dat.imp$treat_lender_pressure
n4b <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
n4 <- ictreg(formula(paste0('list_welfare_pressure ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')
n5b <- ictreg(formula(paste0('list_lender_pressure ~ ', spec1)),
              data = dat.imp, treat = "treat", J = 3, method = 'lm')
n5 <- ictreg(formula(paste0('list_lender_pressure ~ ', spec2)),
             data = dat.imp, treat = "treat", J = 3, method = 'lm')


## table like stargazer

models <- list(p4b, p4, p5b, p5, n4b, n4, n5b, n5)
tab <- cbind.data.frame('mod1' = rep(NA, length(p5$par.treat)*4),
                        'mod2' = rep(NA, length(p5$par.treat)*4),
                        'mod3' = rep(NA, length(p5$par.treat)*4),
                        'mod4' = rep(NA, length(p5$par.treat)*4),
                        'mod5' = rep(NA, length(p5$par.treat)*4),
                        'mod6' = rep(NA, length(p5$par.treat)*4),
                        'mod7' = rep(NA, length(p5$par.treat)*4),
                        'mod8' = rep(NA, length(p5$par.treat)*4))
modinf <- rbind.data.frame('N' = rep(NA, length(models)))
rownames(tab)[seq(1,length(p5$par.treat)*4,2)] <- c(paste0(names(p4$par.treat),".t"), 
                                                    paste0(names(p4$par.control),".c"))

ptmax <- length(p5$par.treat)*2

for (i in 1:length(models)){
  mod <- models[[i]]
  pt <- length(mod$par.treat)*2
  mod$pval <- coeftest(mod)[,4]
  mod$stars <- ifelse(mod$pval<0.01, "***", ifelse(mod$pval<0.05, "**", ifelse(mod$pval<0.1, "*", "")))
  tab[seq(1,pt,2), i] <- paste0(round(mod$par.treat, 3), mod$stars[1:length(mod$par.treat)])
  tab[seq(ptmax+1,ptmax+pt,2), i] <- paste0(round(mod$par.control, 3), mod$stars[(1+length(mod$par.treat)):(length(mod$par.control)*2)])
  tab[seq(2,pt,2), i] <- paste0("(",round(mod$se.treat,2),")")
  tab[seq((ptmax+2),ptmax+pt+1,2), i] <- paste0("(",round(mod$se.control,2),")")
  modinf[i] <- dim(mod$data)[1]
}
colnames(modinf) <- colnames(tab)
xtable(rbind(tab,modinf))






